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Abstract 

We  consider  a  finite  difference  approximation  to  the 
Gauchy  problem  for  a  first-order  hyperbolic  partial  differ- 
ential equation  using  different  mesh  spacings  in  different 
portions  of  the  domain.   By  reformulating  our  problem  as  a 
difference  approximation  to  an  initial-boundary  value 
problem,  we  are  able  to  use  the  theory  of  H.  0.  Kreiss  and 
S.  Osher  to  analyze  the  L„  stability  of  our  scheme. 

In  the  case  of  one  space  dimension,  a  proof  is 
presented  for  the  L  stability  of  the  Lax-Wendroff  scheme 
when  the  interface  condition  is  generated  so  as  to  give 
second  order  accuracy.   Stability  is  also  proven  for  all 
3-point  dissipative  schemes  which  are  matched  along  an 
interface  by  various  conditions. 

For  the  case  of  two  space  dimensions,  stability  is 
proven  with  an  analogous  mesh  refinement  performed  along  a 
coordinate  line.  Along  this  line  data  is  generated  by 
interpolation. 

Numerical  experiments  were  performed  on  the  CDG-66OO. 
We  computed  test  problems  for  the  schemes  discussed  in  our 
theoretical  treatment,  but  also  introduced  an  additional 
feature.   We  used  uneven  time  steps  in  the  different  mesh 
patterns.   Our  numerical  results  indicate  that  our  techniques 
are  an  efficient  method  for  improving  accuracy  in  regions  of 
special  interest,  while  preventing  new  inaccuracies  at  the 
interface. 
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INTRODUCTION. 

There  are  many  situations  in  the  numerical  solution 
of  partial  differential  equations  by  finite  difference 
approximations  where  one  wishes  to  perform  the  calculations 
with  different  mesh  spacings  in  different  portions  of  the 
domain.  This  paper  considers  the  !„  stability  of  a  finite 
difference  solution  to  the  Cauchy  problem 

u,  =Au+Bu    ,  t>0 

t      X      y   ^  — 

(0.1) 

u(x,y,0)  =  *(x,y)  ,  -oo  <  x,y  <  co 

where  the  spatial  mesh  sizes  have  been  refined  in  a 
portion  of  the  domain.   Problems  of  this  type  shall  be 
referred  to  as  mesh  refinement  problems  in  this  paper. 

Chapter  I  contains  a  review  of  the  terminology 
and  basic  theory  of  stability  for  pure  initial-value  and 
initial-boundary  value  problems.   A  specific  mesh 
refinement  problem  is  formulated  for  the  one  dimensional 
case  of  (0.1).   Stability  is  analyzed  by  employing  the 
notion  of  local  normal  modes  as  described  by  Richtmyer 
and  Morton  [11^12]  in  their  exposition  of  the  work  of 
Godunov-Ryabenkii  [1,2].   The  technically  crucial  result 
in  Lemma  1  states  that  the  modes  of  dissipative  equations 
are  separated  by  the  unit  circle. 

In  Chapter  II  the  problem  is  reformulated  as  a 
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difference  approximation  to  an  initial-boundary  value 
problem  for  a  system  of  partial  differential  equations. 
The  bo\xndary  condition  employed  is  derived  from  the 
particular  approximation  used  to  define  the  difference 
equation  at  the  point  of  contact  between  the  two  grid 
patterns.   General  sufficient  conditions  for  the  stability 
of  such  systems  are  due  to  H.  0.  Kreiss  [6]  and  S.  Osher 
[9,10].   In  Chapter  II,  we  review  those  aspects  of  the 
Kreiss-Osher  results  which  we  use  in  this  paper. 

In  Chapter  II,  for  the  case  of  one  space  variable, 
we  show  the  Lp  stability  of  a  generalization  of  the 
Lax-Wendroff  scheme  to  the  mesh  refinement  technique. 
In  Theorem  IV,  we  give  a  proof  of  the  stability  of  all 
3  point  dissipative  schemes  to  the  one  dimensional  case 
of  (O.l)  using  various  interface  conditions.   This  result 
allows  us  to  match  two  schemes  in  a  stable  way  along  an 

interface. 

The  case  of  two  space  variables  is  treated  in 
Chapter  III.   Using  the  results  of  Kreiss  and  Osher,  we 
are  able  to  prove  stability  when  the  mesh  refinement  is 
performed  along  a  line  parallel  to  an  axis,  when  our 
interface  condition  is  generated  by  interpolating  the  data 
along  this  line. 

Numerical  results  obtained  on  the  CDC-6600  are 
presented  in  Chapter  IV.   We  use  these  results  to  compare 
the  accuracy  of  the  various  approximations  discussed. 
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The  numerical  results  bear  out  our  theoretical  conclusions, 
and,  furthermore  yield  evidence  of  stability  for  cases  not 
analytically  treated  (e.g.  variable  coefficients,  refine- 
ments having  comers,  and  refinements  using  different  time 
steps  in  different  portions  of  the  domain). 
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CHAPTER  I.   STABILITY  THEORY. 

a.    Initial-value  problems. 

We  consider  the  numerical  solution  of  the  initial- 
value  problem 

u .  =  A  u       f  -co  <x<oOj   t>0; 

(1.1) 

u(x,0)  =  <t(x)   , 

where  A  is  a  symmetric  matrix,  and  <t>(x)  is  the  prescribed 
initial  value  function  of  the  vector  valued  solution  u(x,t) 
Let  <t)(x)  be  s\iff iciently  smooth  so  that  (l.l)  possesses  a 
solution  operator  which  is  bounded  in  Lp  for  all  finite  t. 
In  such  a  case,  we  say  that  the  problem  (1.1)  is  properly 
posed. 

We  introduce  a  set  of  uniform  mesh  points 

X  =  jh  for  J  =  0, +1, +2, ... 
(1.2)  R=   j(x,t): 

^       t  =  nk  for  n  =  0,  1,  2j  .  ..•' 

where  h  =  Ax  Is  the  spatial  mesh  size  and  k  =  At  denotes 
the  time  step.   Throughout  this  paper  we  assume 

k/h  =  constant. 

Consider  an  explicit  finite  difference  approximation  to 
(1.1)  of  the  form 
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v(x,t+k)  =   Sy^   v(x^t)   =   >   s   v(x+jh,t) 


J 
(1.5) 

v{x,0)      =   u(x,0) 


where  the  s  .  are  constants,  and  (x,t)  e  R. 
J 

When  reference  is  made  to  a  norm  in  this  paper, 
we  shall  always  mean  the  Lp  norm. 


(1.4)  llu(t)|l   =   /  ^~   |u(x,t)|2  ;   t  fixed. 

y  (xTTJeR 

Definition.   The  difference  scheme  (1-5)  approximates  (l.l) 
with  m-th  order  accuracy  if  for  all  sufficiently  smooth 
solutions  u(x,t)  of  (1.1) 

(1.5)  llu(x,t+k)  -S^u(x,t)|!   =   0(h"^+^)  . 

S,     is    said   to  be   consistent   if  m  >  1. 
h  — 

Definition.      The   difference   scheme    (1.2)    is   stable   if 

(1.6)  ||sj|||   <   K  , 

K  some  constant  for  all  n,k  satisfying  n-k  <^  1. 

The  following  theorem  is  well  known  [12]: 
Convergence  of  a  consistent  difference  scheme  to  the 
solution  of  a  properly  posed  initial-value  problem  is 
equivalent  to  the  stability  of  the  difference  scheme. 

Schemes  with  accuracy  of  order  m  can  be  constructed 
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with  the  aid  of  Taylor's  theorem.   For  u(x,t)  a 
solution  to  (1.1)  where  A  is  independent  of  t,  we  can 
express  the  time  derivatives  as  spatial  derivatives  in 
the  Taylor  expansion,  so  that 

(1.7)  u(x,t+k)=  XZ^  (A  ^)^  u(x,t)  +  0(k"^+l)  . 

In  (1.5)  we  expand  v(x+ jh,t)  in  a  Taylor  series  about  x, 
Equating  terms  of  ordSr  m  or  less  in  (1.3)  and  (1.7) 
yields  a  system  of  m  +1  equations  (powers  of  h)  for  the 
m  +1  unknowns,  s„jS,,...,s  ,  the  coefficients  of  S, 
in  (1.3). 

An  alternate  approach  is  to  verify  (1.5)  for   a 
dense  set  of  smooth  data.   One  of  the  advantages  of  this 
method  is  that  one  is  led  to  consider  the  following 
fionction  related  to  (1.3) 

(1.8)  C(|)   =  ZZ  ^•  e 

<]   ^ 

called  the  symbol  (amplification  matrix)  of  the  operator 
S,  .   As  shown  in  [7]?  we  can  study  stability  of  S, 
explicitly  in  terms  of  the  xiniform  boiindedness  of  the 
matrices  u  (^). 

A  well  known  necessary  condition  for  stability 
is  the  von-Neujnann  condition:   the  eigenvalues  of  C  do 
not  exceed  one  in  absolute  value. 

In  the  future  we  shall  refer  to  the  following 
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initial-value   problem  as! 

Exsunple    I. 

u,      =     a  u      ;,  -oo    <  X   <  oo 

(1.9) 

U(X;,0)  =  <J)(X) 

where  a.  is  a  constant  and  u  and.  <t>  are  scalar  functions. 
One  can  verify  that  the  consistency  conditions  for  (I.3) 
give   for  b  =  a  At /Ax 

>   J^  s.  =  b^  ,         i-0,l....,m. 
J      ^ 

If  we  use  a  symmetrically  centered  scheme  of  2N  + 1  points, 
we  are  led  to  take  m  =  2N.   For  the  case  N  =  1,  the 
Lax-Wendroff  scheme   gives  us  second  order  accuracy  using 
as  few  points  as  possible.   For  future  reference  we  note 
that  for  this  scheme 

/.  .  ^x  b^-b  -,    ,2  b^+b 

(1.10)  s_^  =  — 2—  ,        Sq  =  1  -  b   ,    s^  =  — 2— 

where  b  =  aAt/Ax   is  a  constant. 

Using  (1.10)  we  find  as  in  [8]  the  following  equation 

for  the  absolute  value  of  the  amplification  factor  squared, 

(1.11)  |C(^)|2  =  1  -  (b2-b'^)(l-  cos  ^f    . 

Thus  the  von-Ne\imann  condition  (which  is  also  a  sufficient 
condition  for  stability  in  the  scalar  case  [12])  is 
satisfied  iff  Ibi  <  1. 
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b.    A  simple  mesh  refinement  problem. 

Consider  the  Cauchy  problem  of  Example  I.   Suppose 

that  we  attempt  a  numerical  solution  using  two  different 

patterns  of  mesh  spacings.   On  the  right  hand  side  of  the 

origin  the  mesh  length  is  Zix-,  j   on  the  left  hand  side,  we 

take  the  spatial  mesh  size  to  be  Axp.   Setting 

Ax 
(1.12)  p=^ 

we  assume  that,  say,  the  mesh  on  the  right  hand  side  is 
more  refined  than  on  the  left,  i.e. 


0  <  p  <  1  . 


Set 


(1.13)      ^=^    and    ^2  =  2^ 

so  that 

^2  "  P  ^1  • 

Let  w  .  denote  our  difference  approximation  to 

J 

the  solution  u(x,t)  of  Example  I  at  x  =  jAx-, ,  t  =  nAt 

for  j  =  0,1,2,...;  n  =  0,1,2,...  .   Then  for  j=  1,2,5,..., 

2  2 

M  n;no   ,  n+1    /^1"^1n  n    ,  /,  ,  2x  n  ,  /^I'^^lv  n 
(l.l4)a   Wj    =  ( — 2 — )Wj_i  +  (l-b^)w^  +  ( — 2 — )Wj^_i  • 


The  above  formula  also  applies  as  an  approximation 
to  u(j  Ax  ,  n  At)  for  j  negative  if  we  replace  b,  by  b  . 
But,  regardless  of  the  approximation  chosen  for  w^"*"  ,  the 
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difference  equations  have  a  discontinuity  at  x  =  0  in 
the  sense  that  they  approximate  the  differential  equation 
(1.9)  where  the  coefficient  a  =  a(x)  is  a  piecewise 
constant  function  with  a  jump  at  x  =  0.   We  avoid  discuss- 
ing this  by  reducing  our  problem  to  an  initial-boundary 
value  problem  in  the  quarter  plane.   To  this  end,  we 
introduce  a  new  difference  function  v^  which  denotes  an 
approximation  to  u(-j  Ax^,  n  At)  for  j  >^  0.   Then  for 
j  =  1,2,J),  .  .  . ,    and  n  =  0,1,2,..., 


(I.l4)b   v' 


.n+1 


( 


^2+^2 ^  n 


2x  n 


b2-b 


UU    ^    ^l-^2)-i  +  ( 


2   2n  n 


2   /'j_i  ■  ^---2^"j  '  ^-^r-)v 


J+1  • 


We  introduce  the  vector  notation 


(1.15) 


vH 


W  . 
J 

V^ 


;   j  =  0,1,2,...;  n  =  0,l,2,.... 


and  note  that  by  our  definitions 


n  _  n 
^0  =  ^0  • 


Then  (l.l4)a,b  can  be  expressed  as 


(1.16) 


Z^+1  =   Q  Z^ 


where 
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'^'-1 


Q 


(1.17) 


0 


0 

2  2 
P  b^+pb^ 

2 


„-l 


+ 


\ 


1-b- 


0 


0 


2  2 
1-p  bi 


+ 


/  2 

2 

0 


v 

0 
2  2 


,0 


k 


where  we  have  used  (I.I5)  and  where  the  T  ,    the  transla- 


i^n 


^n 


tlon  operators  J  are  such  that  T  Z .  -  ^.- ,  «• 

There  remains  to  prescribe  an  approximation  to 
u(0, t)  at  each  time  step.  In  general  this  would  take 
the  form 


(1.18)   w 


n+1 
0 


fc'      n  n 

y        a.  w.  +  y        P-  V.  ;    s,t   integers 


J  =  0 


0=1 


We  consider  the  case  0  <  p  <  1,  where  we  make  the 
natural  requirement  that  our  difference  scheme  be  of 
second  order  accuracy  at  the  point  x  =  0.   Using  the 
consistency  conditions  derived  from  Taylor's  theorem, 
we  solve  for  aj-j^a-,,  p]_  in 


(1.19) 


n+1 
0 


w 


OqW^  +  a-j_w^  +  p^v^  +  O(Ax^) 


an  unevenly  spaced  scheme.   A  straightforward  calculation 
yields  the  following  set  of  equations; 
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Uq  +  a-,   +  Pn   =   1 


2 


(1.20)  a-^p  -  p-j^   =  b 

2  2 

Solving  and  substituting  into  (I.19)  we  have  up  to 
second  order  accuracy 

n+1     V^2_lPl  n  ^  ..    ^2w,  ,^  ,  n 

^0   =   (1  +  p)  ^1  +  (1  -  -p-)  (1+^2)^0 

(1.21) 

VV^  n  . 
+   p(p+l)  ^1  ^   ^  =  0^1^2,. 

We  can  express  (1.21)  in  the  following  vector  form 
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absolute  value  then  we  expect  that  one  could  construct 
an  unstable  mode  ( ]  A]  >  l)  by  using 


,   _  /  "l^i   +   ^2^^2 


(1.26) 

J      \         H-v, 

which  could  be  chosen  to  satisfy  the  boiindary  conditions 
(1.22).   That  this  cannot  be  the  case  is  due  to  the  fact 
that  the  ^-^y^n    ("^i^"^?^  ^^^   separated  by  the  unit  circle 
so  that  always,  say,  JM--,]  <  1  and  \l^A    >   1  for  |  a|  >  1. 


d.    Separation  oC  the  modes. 

Consider  a  consistent  three  point  difference 
approximation  to  (1.9)  for  n  =  0,1,2,...;  j  =  0,+l,+2, . , 

(1.27)         wj   =  a_^w^._^  +  a^w^.  +  a^w^^^ 

a  At 
where  the  a.  are  constants.   Let  b  =  ^—  be  a  constant, 

We  check  the  von-Ne\imann  condition  by  substituting  into 

(1.27)  a  solution  of  the  form 

(1.28)  w^  =  c  tx^A^  ,     c,  constant, 

where  m-  =  e"""  ,  ^  real,  |  ^|  <  tt. 

We  get  a  quadratic  equation  in  I-l  for  the  modes; 

(1.29)  a^ix^  +  (a^-  A)li.  +  a_^  =  0  . 
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The  difference  scheme  (1.27)  is  stable  iff  ]  p.]  =  i  implies 
I  a|  <  1. 

Lemma  1   (Separation  of  roots) . 

Let  the  difference  scheme  (1.27)  be  stable^  so  that 

(1.30)   |a_^  e"-^^  +  a^  +  a^  e^^|  ^1,   |4|  <  tt  ,  ^  real. 

Then  the  \i   roots  of  (1.29)  are  separated  for  any  A, 

I a|  >  1,  which  means  that  one  of  the  roots  of  (1.29)  lies 

inside  the  unit  circle,  and  the  other  outside.   (Note, 

more  general  results  of  this  form  can  be  found  in  [6]  and 

[9].) 

Proof.   Let  V^-^,    V-^   be  the  roots  of  (1.29),  then 

a^  -A 

^1^2  "  ^  l/^l     ^^^      \^^   ^  "^^  =  -^ . 

"^1 

Taking  the  limit  as  A  tends  to  infinity,  we  observe 
that  the  product  remains  constant,  but  the  sum  blows  up. 
Thus,  we  conclude  that  one  root  tends  to  infinity  and  the 
other  to  zero.   That  is,  the  unit  circle  separates  the 
roots  as  A  — >  oo  .   That  this  is  true  for  any  A„  such  that 
|a  1  >  1  follows  from  the  fact  that  the  roots  M--,(A), 
i  =  1,2,  depend  continuously  on  A.   Indeed,  as  A  tends 
from  infinity  to  A  through  values  of  absolute  value 
greater  than  one,  the  values  M--,  ( A) ,  l-i-p(A)  cannot  touch  the 
unit  circle,  for  then  the  stability  criterion  would  imply 
that 
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absolute  value  then  we  expect  that  one  could  construct 
an  linstable  mode  ( j  a|  >  l)  by  using 

/c  M-^'   +  c  M-J 
(1.26)         i_   =  M  ^    _   ^  ^ 

^        \  dv^ 

•which  could  be  chosen  to  satisfy  the  boundary  conditions 
(1.22).   That  this  cannot  be  the  case  is  due  to  the  fact 
that  the  M--, ,  l-i-p  (v-,jV  )  are  separated  by  the  unit  circle 
so  that  always,  say,  \^-,\    ^  ^   and  1m-„1  >  1  for  |  a|  >  1. 


d.    Separation  of  the  modes. 

Consider  a  consistent  three  point  difference 
approximation  to  (1.9)  for  n  =  0,1,2,...;  j  =  0,+l,+2, . , 

/  n  r^r^\  n+1       n    ,    n  ,    n 

(1.27)         Wj.   =  ^-i^j-i  +  ^o^j  +  ^l^j+1 

a  At 
where  the  a.  are  constants.   Let  b  =   -r —  be  a  constant. 
1  Ax 

We  check  the  von-Neumann  condition  by  substituting  into 

(1.27)  a  solution  of  the  form 

(1.28)  w.  =  c  H-^  A  ,      c,  constant, 

where  I-l  =  e  ,    i   real,  j  ^1  jl  "Ji". 

We  get  a  quadratic  equation  in  l-i-  for  the  modes; 

(1.29)  a-|_|x^  +  (a^-  A)|i  +  a_^  =  0  . 
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The  difference  scheme  (1.27)  is  stable  iff  |  ii.|  =  l  implies 
I  a|  <  1. 

Lemma  1   (Separation  of  roots) . 

Let  the  difference  scheme  (1.27)  be  stable,  so  that 

(1.30)   |a_-j_  e"-^^  +  Bq  +  a-j^  e^^|  ^1,   Ul  1  "^  ,  ^  real. 

Then  the  14.  roots  of  (1.29)  are  separated  for  any  A, 

I  A|  >  1,  which  means  that  one  of  the  roots  of  (1.29)  lies 

inside  the  unit  circle,  and  the  other  outside.   (Note, 

more  general  results  of  this  form  can  be  found  in  [6]  and 

[9].) 

Proof.   Let  l^-^,    M-^  be  the  roots  of  (1.29),  then 

^0  "^ 

^1^2  "  ^_l/^l      ^^^  M-n  +  M-o  =  . 

"^1 

Taking  the  limit  as  A  tends  to  infinity,  we  observe 
that  the  product  remains  constant,  but  the  sum  blows  up. 
"Thus,  we  conclude  that  one  root  tends  to  infinity  and  the 
other  to  zero.   That  is,  the  unit  circle  separates  the 
roots  as  A  — >  oo  .   That  this  is  true  for  any  A„  such  that 
|a  I  >  1  follows  from  the  fact  that  the  roots  |x  (A), 
i  =  1,2,  depend  continuously  on  A.   Indeed,  as  A  tends 
from  infinity  to  A  through  values  of  absolute  value 
greater  than  one,  the  values  M--,  ( A) ,  M-p(A)  cannot  touch  the 
unit  circle,  for  then  the  stability  criterion  would  imply 
that 
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I  A]  <  1  . 

Thus  M--,  and  M-  remain  separated  by  the  imit  circle  for 
all  A,  I  a|  >  1.  Q.E.D. 

By  consistency,  a_-,  +  a^  +  a-,  =  1.   So  M-  =  1  implies 
A  =  1.   However,  the  following  is  easily  proven  by  a 
similar  arg\ament  as  in  Lemma  1. 

Corollary  1.   If  for  0  <  \^\    _^  tt  equation  (I.50)  is 
satisfied  with  strict  inequality  then  the  unit  circle 
separates  the  two  M-  roots  of  (1.29)  for  j  AJ  _>  1,  A  7^  1. 
We  see  from  (l.ll)  that  the  Lax-Wendroff  scheme 
satisfies  the  hypothesis  of  the  above  corollary,  and  thus 
all  Lp  solutions  of  the  form  (1.23),  (1-24)  to  (I.I6), 
(1.17)  must  have  the  simple  form 

/  C  ^L^"  \ 

(1.24)'        in-  =       i  K   J  =0,1,2,...;  c,d,  const. 

J     \  d  v^  / 

A  necessary  condition  for  stability  of  our  mesh  refinement 
problem  of  Section  b.  is  that  (I.23),  (1.24)'  does  not 
satisfy  (1.22)  for  | a]  >  1. 

It  can  be  shown  that  the  interface  condition 
according  to  (1.22)  rules  out  the  possibility  of  having 
1  a]  >  1  for  our  solution  (1.24)'.   No  verification  is 
made  here,  because  this  result  will  be  a  consequence  of 
our  proof  of  sufficiency  conditions  for  stability  in 
Theorem  I. 
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CHAPTER  II.   MESH  REFINEMENTS  AS  INITIAL- BOUNDARY  VALUE 
PROBLEMS . 

a.    The  theory  of  Krelss  and  Osher. 

In  this  section^  we  present  the  theory  of  Krelss  [6] 
and  Osher  [9^10]  for  initial-hoimdary  value  problems 
specialized  to  a  simple  case  which  suffices  for  our  purposes, 
Consider  a  first-order  system  of  partial  differential 
equations 


Su 


Su 


/  „    -,  X       du         .    du     ,    T-,   du         ,  _         _ 

^     ^^      ^  ^        '5x  ^  '    t>0,    0^x<oo    ,    -00    <y<oo, 

where  A  and  B  are  constant  diagonal  matrices .   A  has  the  form 


r  a- 


(2.2)    A 


a 


2 


^ 


a 


K 


'K+1 


a 


a-,  <  0,  . .  .  ,a„  <  0 


^K+1  ^  0,  .  .  .  jBj^  >0. 


M 

We  consider  an  initial-boundary  value  problem  in  the 
quarter  space  x  >  0^  t  >  0. 


(2,3)   u(x,y,0)  =  Mx^y)  ,       0<x<oo,   -oD<y<oo. 


(2.4)a)    u^(0,y,t)  =      S^^  u^^(0,y,t)  , 

where  S    is  a  constant  rectangular  matrix^  and  where 
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(2.4)  b) 


u^  =  (u(l)(x,y,t),...,u(^)(x,y,t)) 


II        I         (K+1),  ,N  (M)/  ,SV 

u   =  (u^    ^(x,y,t),  .,  .,u^  Mx,y,t)) 


We  approximate  (2.1)- (2.4)  by  the  system  of  difference 
equations 

for  n  =  0,1,2,...;  i   =  s,s+l,...;  j  =   ...  -1,0,1,... 
with  boundary  conditions 

n      r+s   ^2        i^  ■ 

for  i  =  0,-1, .  .  .,-s+i;  j  =  ...  -1,0,1,...;   n  =  0,1,2,... 
where  q,r,r-, ,  rp,  s,  t-,  ,  tp  are  all  non-negative  integers,  and 

(2.6)  ^\^^    =  (v'^^j(l),...,v'^^^.(M)) 

is  an  M  component  vector  approximation  to  u(  ^Zix,  jAy,nAt ) . 
The  C   Q  are  constant  diagonal  (M  x  m)  matrices  and  the 
b.     are  constant  (M  x  m)  matrices. 

We  assume  that  a  time  step  k  =  At  >  0  and  mesh 
sizes  h-,  =  Ax  >  0  and  hp  =  Ay  >  0  have  been  introduced, 
so  that  b-,  =  k/h-,  and  bp  =  k/hp  are  constants. 

The  following  three  conditions  are  natural 
assumptions  for  the  stability  of  (2.5)1 
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1)  The  difference  approximation  (2.5)  with  its 
boundary  conditions  is  consistent  with  the  initial- 
boundary  value  problem  for  the  differential  equation  (2.1), 

2)  The  difference  approximation  defined  by  (2.5)a 
is  Lp  stable  in  the  pure  Cauchy  case. 

3)  The  difference  operator  C  =  ">   C   ^  rp-^+a  rjnJ+P 
defined  by  (2.5)  is  dissipative  [5]  in  the  following  sense! 
the  eigenvalues  A  =  A(^)  of  the  amplification  matrix 

G  =  ^:  C   ^  e^((^+°')^l+(j+P)^2)^  the  Fourier  transform 
of  C,    are  strictly  less  than  one  for  c,   t,   real  and 
0  <  Al+i^  ITT. 

A   sufficient  condition  for  stability  is  found  in! 

Theorem  (Kreiss-Osher) . 

Suppose  that  the  assumptions  I-5  are  valid. 
Suppose  further  that  (2.5)a,b  have  no  non-trivial 
generalized  solutions  of  the  form 

(2.7)  v^^.  =  A^I.yj  ,     IyoI  =  1  > 
for  which  either  a)  or  b)  or  c)  holds^  where 

a)  1  a]  >^  1^  A  /  1  and  y„  7/  1  and  >    |  0^  J   <  oo 

^         £=0         ^ 

00       2 

b)  I  a|  =1,  A  ^  1,  y  :=  1    and  ^^   \^A       <  00 

(2.8)  _  '-° 

c)  A  =  y  =  1  and  1*^.|  _<  constant, 

i-^I|  1  li        II 

\^p\    Jl  constant  ]  M-|   where  |  M.|  <  1. 
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Then  the  approximation  (2.5)  is  stable.   (See  (2.4)b 


for  the  definition  of  ^., 


) 


Note  that  the  above  theorem  holds  for  the 
one-dimensional  case  B  =  0  with  the  obvious  simplification 
Yq  s  1,  [6,9]. 


b.     Stability  of  a  mesh  refinement  problem. 

We  return  to  the  one-dimensional  mesh  refinement 
problem  of  Chapter  I,b.   Equations  (l.l'^  )- (1. 1?)  remain 
with  only  one  minor  change.   As  before,  w.  denotes  an 
approximation  to  u(j  Ax-.,  n  At)  for  j  =  0,1,2,...; 
n=  0,1,2,...,   but  now  in  order  to  conform  with  the 
formulation  of  Kreiss  and  Osher,  we  found  it  necessary 
to  change  our  notation  so  that  v^  denotes  an  approximation 
to  \x{-{j-l)I^^,    n  At)  for  j=  0,1,2,...;  n=  0,1,2,...  . 
(See  figure  below.) 


t  t 

At 

4  4 


Ax, 


W- 


w. 


^iSx  -*■ 


"0=^1 


V, 


Figure  1 
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From  our  definitions  of  v  .  and  w  .   Vie  have 


2.9) 


w 


n 

0 


n 


V- 


n  =  0,1,2, 


In  general,  a  boundary  condition  of  the  type  (2.5)b 
represents  an  extrapolation  along  a  given  time  step.   What 


n 


is  missing  in  our  case  is  a  formula  for  v„.   We  obtain  this 
by  second-order  extrapolation  along  the  x-axisi 


(2.10)   Vq  =   a^wj  +  a^v^   +   a^v^  +  O(Ax^)  . 


n 


n   n 


Expanding  w, ,  v-,  ,  v„  about  the  interface  and  using  (2.10) 
to  achieve  second  order  extrapolation,  we  obtain  the 


following  equations  for  the  a. , 


a-,  +  Qp 


(2.11)<  (l-p)a.    +  a 


1 


(1-p)' 


V.   2 


a-j^  +  2  ct2 


+ 


a 


5 


+    2  a. 


3 


+    2  a- 


=  1 


=  0  where  p 


=  0 


Ax 


1 


Ax, 


Solving  (2.11)  we  get 


(2.12)   a^  = 


,  g  _  2(p-l)       _  !-_£ 

pTT+pT  '  "^2  -   p   ^  "3  ~  1  +  p  • 


We  take  as  a  defining  approximation  for  n  =  0,1,2,..., 
/„  -i-jN    n      2     n  ,  2(p-l)   n  ,  /l-p>,   n 
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th  at  V' 


Using  the  above  equation  for  v^  in  (I.I6),  we  find 
n+1 


1 


is  then  determined  as 


'T' 


(b^+b^) 


-I ^wj  +  (1-b  /p)(l  +  b  )v5 

(1+p)    -*-        ^        ^  -^ 


(2.14) 


+ 


VV^P^^n 


(1+p) 


2  ' 


n  =  0,lj  2,  .  . .  . 


When  one  looks  at  the  definitions  for  v.  and  w.  in 

J      J 

Chapter  l,h   and  here,  the  above  expression  for  v-j^   is 
seen  to  agree  exactly  with  equation  (1.21),  which  is  a 


direct  derivation  of  v' 


n+1 
1 


as  a  linear  combination  of 


w 


n 


n 


V-, 


n 


Vr 


retaining  second-order  accuracy  at  x  =  0. 


1'  "1^  "2 
Thus  the  formulation  of  this  section  is  equivalent  to  that 

of  Chapter  I,b. 

In  order  to  apply  Kreiss'  theorem,  we  need  only  give 

a  vector  formulation  of  (2.9)  and  (2.15)  which  will  serve 

as  our  boundary  condition  for  (I.I6),  (I.I7).   Namely, 


(2.15)  z", 


0 


f.l^ 


n 


0 


V°i 


a 


\ 

«i 

+ 

i 

2) 

n 

K 

0 


0 


0 


^  '  ^2 


a. 


5/ 


n 

V^2  ; 


n 


0,1,2, 


where  the  a.  are  defined  in  (2.12).   We  note,  that  because 
our  difference  equations  were  obtained  by  reflection  of 
the  left  hand  side,  (I.I6),  (l.l?)  can  be  shown  to  be 
consistent  with  the  system 
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(2.16) 


CI 


a 
.0 


0  ^ 
-ap. 


w\ 


Note  that  the  above  diagonal  elements  are  always  of 
opposite  sign.   This  will  he  very  pertinent  when  we 
wish  to  verify  condition  c)  of  Kreiss'  theorem. 
We  shall  refer  to  the  system  of  difference  equations 
(l.l6)j  (1.17)  with  boundary  condition  (2.I5)  as 
System  I. 

Theorem  1.   System  I  is  a  stable  approximation  for  the 
Cauchy  problem  of  Example  I  for  jb-.  |  <  1. 

Proof.   Let  a  <  0.   A  similar  analysis  can  be  employed 
for  the  case  a  >  0. 

As  seen  in  Chapter  I,c,d^  we  find  all  solutions 
to  (1.16),  (1.17)  of  the  form 


,n    -,n  -r 


00 


(2.17)a)   z^j  =  A^^  £.  ;  ZZ\1^\      ^"^^      A^l,  1A|>-1 


■J   j=o   J 


by  taking 


(2.17)t)) 


(^1  4  \ 

*  .    = 

~J 

P2  i  , 

p.   constants J 


V 


where  r   ,x     are  determined  respectively  from 
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2  2 

a)  A  =  (^^)^  +  (l-b^)  +  (-i_l)  T^ 

(2.18)  22'  22 

P  b,+pb   .  p  b  -pb 

b)  A  .  ( V^)l^  +  (l-P^b^)  +  ( 1^)    T^ 

In  Chapter  I,C;,d  we  showed  that  T  .  was  of  the  form  (2.I7) 
by  showing  that  the  modes  of  the  Lax-Wendroff  scheme  were 
separated  by  the  ■unit  circle  for  |  Aj  >^  1,  A  :j/  1.   By 
Kreiss'  theorem^  System  I  is  stable  if  it  has  no  solution 
of  the  type  (2.I7)  nor  of  the  type  (2.8c),  that  is_, 

(2.19)  Z.  =  "*.  where  I*  .  I  <  constant  and  ' 

I"*.  I  <  const. Im-H,  I  uI  <  1. 

For  the  case  A  =  1,  we  can  solve  (2.l8)  directly  for 
T.   ,  i  =  1,2.   Observe  that  in  this  case,  since  a  <  0, 

I T   1  >  1   and   1 Tp  I  =  1    . 

Referring  to  (2.19),  we  can  exclude  the  above  modes  from 

■^..   Thus  for  the  case  A  =  1  the  form  of  the  *".  is  still 
-J  -J 

(2.17).   Because  a  <  0,  checking  (2.8c),  we  see  that  an 
unstable  mode  must  be  of  the  form  (2,17)  with 


We  now  prove  stability  by  showing  that  all 
non-trivial  solutions  of  system  I  of  the  form  (2.I7)  must 
have   T  =  1. 
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Substitution  of  (2.I7)  into  the  boundary 
condition  (2.I5)  yields 

^)   ^1  ~  P2'^2 
(2.20) 

b)   ^2  =  o'lPl^l  +  a2^2'^2  +  «3p2'^2 

where  the  t^  are  defined  by  (2.18)  and  the  a.  are  defined 
by  (2.12).   ¥e  are  looking  for  a  non-trivial  solution  for 
P^,p2-   Clearly  by  (2.20)a   p^  7^  0  for  a  non-trivial 
solution.   Using  the  first  equation  to  eliminate  p.,  in 
the  second,  and  using  the  values  for  the  a.,  we  find, 
using  that  p  ^  0,    that 

(2.21)         p(l-p)T^  +  2T^T2  +  2(p2-l)T2  -  p(l+p)=  0. 

Equating  the  expressions  for  A  from  (2.I8),  yields 
upon  clearing  fractions 


2b^(l-p2)T^T2  +  T^[p2b^+pb^+(p2b^_pb^)T2_(b^+b^)T^T2] 


(2.22) 


-  (b^-  \)t^      =      0 


Solve  for  ^-j^x  and  x  from  (2.21)  and  substitute  into  the 
above  to  get 
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(2.23)  [P(P+1)  +2(1-p^)t2  +  p(p-1)t2]  ^2hl{l-v^)r^   +   b^p^ 
+  b^p  +  (b2p2_b^p)T2-(b2+b^)[PiMll+(i_p2)^^ 
^2l^rl^]      .2    (b2.b^)x2   .   0. 

The  above  expression  is  a  polynomial  of  degree  4  in  t 
with  coefficients  in  terms  of  p  and  b, .   Let  us  denote 
the  left  hand  side  of  (2.25)  by 

(2.24)  F(p,b,T2)  =  Aq  +  A-^Tg  +  A^Tg  +  A^t|  +  A^^t^  . 

We  can  determine  the  coefficients  A.  by  a  straightforward 

inspection.   We  find  that 

b^-b 
Aq  =  (^^)  P^(P^-I)  .    A-L  =  2p2(i_p2)(b2_^^)  ^ 

(2.25)  A^  =  3(b2-b^)(p2-l)    ,    A^  =  2p2(b^-b^)(l-p2)  , 

A,  =  P!i^  (b^-b,)  . 
Thus  (2.24)  becomes 

F(p,b^,T2)  =  (-2lA)p2(p2_i)(l_4^^^^2_4^5^^^) 

(2.26)  2 

=  (-^)p^(p^-i)(Vi)^  =  0  • 

For  a  reasonable  problem  we  ass-ume  that  p  is  neither 
0  nor  1.   Thus  J  since  we  have  assumed 


(2.27)  0  ^  Pi  I  -^  1  . 


1 
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we  conclude  from  (2.26)  that 


(2.28)  T   =  1  . 

Q.E.D, 


Remarks. 

l)   For  the  case  a  >  0,  we  need  to  alter  the 
above  argioment  in  two  places.   For  A  =  i^  now  t  -=   1  and 
T    >  1.   Condition  c)  of  Kreiss'  theorem  dictates  that 

2+ 

*_.  is  again  of  the  form  (2.17)b.   The  derivation  of  (2.26) 
J 

can  now  be  carried  out  as  in  our  theorem.   But  now,  an 
unstable  mode  may  have  t^  =  1  if  A  =  1.   However,  if 
T  =  1  in  (2.21),  we  see  that  t  =  1.   Since  an  iinstable 
mode  for  this  case  must  be  of  the  form 


we  have  proven  stability  for  the  case  a  >  0. 

2)  The  above  proofs  are  valid  for  the  case  p  >  1, 
since  we  only  used  the  fact  that  p  was  not  0  or  1.   Thus 
System  I  is  stable  regardless  of  which  side  is  refined 
with  respect  to  the  other. 

3)  A  separate  verification  of  Kreiss'  conditions 
shows  the  stability  of  System  1  when  lb-,|  =  1. 
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c.    Boundary  conditions  for  general  mesh  refinement 
problems . 

The  boundary  condition  employed  in  the  previous 
section  for  System  I  was  determined  by  the  requirement 
of  second-order  accuracy  of  the  Lax-Wendroff  scheme  at 
the  interface  of  the  grid  patterns.   The  algebra  of  the 
stability  proof  was  somewhat  involved.   In  this  section, 
we  present  different  interface  conditions  which  have  the 
advantage  that  stability  of  the  scheme  is  easily  demon- 
strated.  An  additional  advantage  is  that  these  methods 
are  extendible  to  higher  dimensions. 

We  make  use  of  the  following  formula  for  the 
linear  extrapolation  of  a  sufficiently  smooth  function 
v(x). 


Lemma  2.   Let  v.  =  v(j  Ax)  then 
J 

(2.29)     V   =  ^   (-1)J+1  {"})    V   +  O(Ax^)  . 
<J  j=l  ^         ^ 

Proof .   Solve  for  an  approximation  to  v^  with  order  of 


accuracy  n-1  in 


n 


(2.30)  v^  .  ^=J"J- 

Then  (2.29)  follows  directly  by  Taylor's  theorem.  (Or, 
see  a  treatment  of  divided  differences  e.g.  Isaacson 
and  Keller  [5].) 
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For  the  case  n  =  3  then 


(2.31) 


V 


0 


3v-j_  -  3V2  +  v^  +  O(Ax^) 


We  now  alter  System  I  so  that 


(2.32) 


V 


n 
0 


w 


n 
0 


n 


0  T  ? 


as  in  Chapter  I^b.   We  then  use  (2.31)  to  define  Vq  and 
w^  for  n  =  0,lj2,...  .   The  resulting  mesh  refinement 
approximation  scheme  will  be  referred  to  as  System  II. 


Theorem  II.   System  II  is  a  stable  approximation  to 
Example  I  when  a  <  0.   If  a  >  0  then  replacing  (2.31)  by 


(2.33) 


n   V  n   ^  n  ,  ,  n 
Vq  =  3w^  -  3^2  +  w^  , 


n 


Ofl,!t}...f 


results  in  a  stable  mesh  refinement  scheme   for  |b,  |  <_  1 


1 


n+1 


Remark.   If  the  resulting  approximation  for  Vq   does 
not  conform  with  the  direction  of  the  propagation  of 
influence,  the  sufficient  conditions  for  stability  are 
not  satisfied.   That  one  has  instabilities  in  such  cases 
will  be  shown  in  numerical  examples  in  Chapter  IV. 


Proof.   In  matrix  language,  our  boimdary  condition  (2.31) 

fo   -A 


has  the  form 
(2.34) 


/  n^ 
n 


V 


vO/ 


(  n\ 

1 
n 


+ 


0  -3 


"2 
n 

v^2/ 


+ 


'0  1 


0   1 


f   n\ 

"3 
n 

V  5/ 


V 
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Systems  I  and  II  have  the  same  underlying  difference 
equations,  and  differ  only  in  their  interface  conditions. 
Thus,  it  suffices  to  check  solutions  of  the  form  (2.I7)  to 

(2.34)  in  applying  Kreiss'  theorem.   Substitution  of  (2.I7) 
into  the  above  boundary  condition  gives 

(2.35)  Pi  =  P2^^''2  "  ^^^2  +  ^2^  =  ^2  • 
For  a  non-trivial  solution  Pp  ^  0  thus, 

(2.36)  (t^-I)^  =0   or   T^  =  1  . 

Checking  condition  c)  of  Kreiss'  theorem,  we  find 
for  a  <  0  that  t^  =  1  is  an  inadmissible  mode,  which  proves 
stability.   If  however  a  >  0  in  the  above,  then  x  =  1 
would  be  an  admissible  mode.   Then  A  =  t  =  1  with 
p  ^   ^     J  Q>     shows  that  the  sufficient  conditions  for 
stability  are  not  satisfied. 

For  the  case  a  >  0,  we  would  use  (2.35)  "to  obtain 
a  stable  mesh  refinement  problem.  Formulating  (2.35)  in 
matrix  language  and  checking  with  a  solution  of  the  form 
(2.17)  gives  us 

(2.37)  Pi  =  ^1(3^1  -3t^  +  tJ)  =  ^2   • 
For  a  non-trivial  solution 

(t^  -  1)^  =  0  or   T-j^  =  1. 
This  proves  stability  since  a  >  0.  Q.E.D. 
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It  is  easy  to  see  that  boundary  conditions  of  the 
form  (2.29)  will  always  give  us  stable  approximations  if 
the  extrapolation  is  in  terms  of  points  which  are  in  the 
domain  of  dependence  of  v^   .   The  stability  of  such 
boundary  conditions  has  been  discussed  by  H.  0.  Kreiss 
in  [4,6]. 

Next  we  turn  to  a  new  type  of  interface  condition, 
suggested  by  H.  0.  Kriess,  which  has  the  advantage  of  being 
stable  regardless  of  the  sign  of  a,  and  is  also  suitable 
for  the  case  of  higher  dimensions.   However,  we  are 
restricted  to  the  case  where 

P  ~  S5c^  "  M 
where  M  is  an  integer.   Assuming  a  definition  for  v.  and 

J 

w.  as  in  System  1,  we  see  that  v„  coincides  with  w.,  for 
J  0  M 

n  =  0,1,...  .   Thus  our  boundary  condition  in  matrix  form 
becomes 

'01' 

n  -  0,1,2,  .  . 

V  ^0  '  ^°  ° 

We  shall  refer  to  the  difference  approximation  (to  Example  l) 
formed  by  using  (2.58)  in  conjiinction  with  the  difference 
equations  (I.I6),  (I.I7)  as  System  III. 

Theorem  III.   System  III  is  a  stable  approximation  to 
Example  I  for  [b-j^l  <  1. 
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Proof.   Again,  because  of  the  separation  of  the  modes j, 
to  check  stability  by  Kreiss'  theorem,  we  try  a  solution 
of  the  form  (2.I7).  'Substitution  into  the  boundary 
condition  (2.38)  gives 


a) 

Pi   =    P2'^2 

(2.39) 

^) 

^2   =    ^1^1 

A   non-trivial  solution  exists  only  if  P-,  and  p„  are 
different  from  zero.   Eliminating  the  p. 's  we  find 

(2.40)         ^2  ""l  =   ^  • 

Since  we  look  for  a  bounded  solution,  the  above  implies 
that 


Regardless  of  the  sign  of  a,    this  assures  stability, 
since  by  condition  c)  of  Kreiss'  theorem  one  of  the  modes 
must  be  strictly  less  than  1  in  absolute  value. 
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d.    Matching  two  difference  schemes  along  an  Interface. 

In  this  section^  we  consider  the  stability  of  an 

approximation  to  Example  I  formed  by  using  two  different 

5-polnt  stable  difference  schemes.   Let  w^  and  v''^  denote 

J      J 

difference  approximations  to  u(j  Ax, ^  n  At)  and 
u(-(j-l)Ax„j  n  At  )  respectively.   Let  S-,  and  Sp  be  two 
difference  operators  so  that 


(2.41) 


a) 


b) 


w  . 
0 


v^+1 


S,w.  =  a  -,w.  -,  +  a„w  .  +  a-,w^,-, 
1  J     -1  J-1    0  0    1  j+1 


^2^5  =  d_iv5_-L  +  d^v^  -f  d^v5^3_ 


n=0,l,2. 


where  a.  and  d.  are  constants.   Assume^  as  In  System  III, 

that 

Ax. 


P  = 


Ax. 


1    1 
M 


where  M  Is  an  Integer.   For  Z.  defined  by  (I.15)j  we  write 

J 

our  difference  equations  in  matrix  form  as 

fs^        0  \ 


(2.42)a 


,n+l 


1 
0    S 


2j 


Z^  ,   j=  l,2,3....;n=0,l,2,  . 


At  the  Interface  of  the  two  meshes  we  employ  the  boundary 
condition  of  System  III, 


(2.42)b 


Z^  - 


fo      l\ 


0    0 


Theorem  IV.   If  (2.42)a   is  consistent  with  the  system 
(2.16)  and  if  for  0  <  \^\    ^  tt   , 


-J>> 


(2.45) 


a)      la_-^  e  ^^  +  a^  +  a-j_  e^  I   <  1  , 


b)      |d_^  e  ^^  +  dQ  +  d^  e^^l   <  1 


then  the  difference  approximation  (2.42)  is  stable. 

Proof.   To  check  stability  by  Kreiss'  theorem,  we  find 
all  solutions  to  (2.42)  of  the  form 


(2.44)         Z^   =  A^ 

J 


i\ 


ri 


^  ,  ,    j  =0,1,2,...;  n=  0,1,2, 
i^H2  '^2 

where  |  AJ  _>  1  and 

a)  A  =  a_j_  /  T^  +  aQ  +  a^T^ 

(2.45) 

b)  A  =  d_^  /  T^  +  dQ  +  d^T^ 

Because  (2.43)  holds.  Corollary  1  implies  that  for 

aI  >  1  ,        A  7^  1, 


the  unit  circle  separates  the  t  and  x  roots  of  (2.45). 
Thus,  as  in  Theorem  III,  when  Z^  is  substituted  into  our 

J 

boundary  condition  we  have 


hil  =  1^21  =   1  • 

We  need  only  check  condition  c)  of  Kreiss'  theorem  to 
complete  our  stability  proof.  Our  boundary  condition 
would  imply  stability  as  before,  if  we  knew  that  the 
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modes  were  separated  for  the  case  A  =  1. 

Case  A  =  1. 

We  give  the  proof  of  separation  of  the  modes  for 

the  case  a  <  0.   A  similar  proof  holds  for  the  case  a  >  0. 

Consistency  means  accuracy  of  at  least  first  order. 

Taylor's  theorem  gives  us  for  the  difference  scheme  (2.4l)a 

a  -,  +  a^^  +  a-.   =  1 

(2.46) 

-  a_-,  +  a-,       =  b,   =  a  At/  Ax-, 

Using  (2.43)a  and  the  above  we  solve  directly  to  see 
that  in  any  case  of  interest 

(2.47)         |a_^/a^|   >  1   . 

The  inequality  expressed  in  (2.47)  assures  us  that 
one  T  root  is  greater  than  1  in  absolute  value.   This 
mode  can  be  excluded  and  thus  the  t  modes  are  separated. 
We  finally  note  that  the  x  =  l  solution  to  (2.45)b 
(with  A  =  l)  is  to  be  excluded  since  a  <  0.   Thus  the 
modfes  are  separated  and,  as  indicated  above,  our 
boiondary  condition  proves  stability  for  this  case.  Q.E.D. 

The  extrapolation  formula  (2.29)^  when  used  so  as 
to  be  consistent  with  the  direction  of  propagation  will 
also  yield  a  stable  matching  of  the  difference  schemes 
of  Theorem  IV. 
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e .    Matrix  case. 

Without  loss  of  generality,  we  let  A  in  (1.1)  have 
diagonal  form  as  in  (2.2).   To  form  a  stable  mesh  refinement 
to  (1.1 ) J  we  could  employ  a  3-point  dissipative  approximation 
in  each  of  the  M  dimensions.   As  in  the  scalar  case,  we 
would  reformulate  the  problem  as  a  system  of  difference 
equations  consistent  with  an  initial-boundary  value  problem. 
The  resulting  difference  equations  for  the  mesh  refinement 
problem  are  in  diagonal  form.   The  local  normal  mode 
analysis  of  the  vector  case  reduces  to  verifying  stability 
in  each  of  the  M  dimensions  in  the  exact  manner  of  the 
scalar  case.   Using  Kreiss '  theorem,  one  easily  verifies 
that  Theorems  I-IV  hold  for  their  matrix  analogs. 
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CHAPTER  III.   EQUATIONS  WITH  TWO  SPACE  DIMENSIONS. 

We  consider  the  analogous  stability  problem  for 

(3.1)         u^  =  a  u^  +  t)  u  ,    a^b  constants^ 

in  the  plane  -00  <  x,y  <  00  ,  t  _>  0.   We  approximate 
the  above  equation  for  the  Cauchy  problem  by  using  a 
difference  approximation  with  a  different  mesh  pattern 
in  each  half  space.   Let  the  plane  x  =  0  be  the  interface 
between  the  two  mesh  patterns.   In  the  right  half -plane 
we  take  mesh  spacings  in  the  x  and  y  directions  to  be 
Ax   and  Ay^  respectively.   In  the  left  half  plane  we 
denote  these  respective  mesh  spacings  by  Axp  and  Ay„. 
We  restrict  our  attention  to  the  case  where 

^^1   1       ,  ^^1   1 

Pi  =  A5^  =  N      ^^^      P2  =  A^  =  M 

where  N  and  M  are  positive  integers. 

Let  W  .  1  and  vT  ,  denote  approximations  to 
■J .  k      J ,  k         ^^ 

u(j  Ax-^,  k  Ay-^,  n  At)  and  u(-(j-l)Ax2,  k  Ay^,  n  At) 

for  n  =  0,1,2,  ...;   j  =  0,1,2,  ...  ;   k  =  0,  +  1,  +  2,  .  .  .  . 

(See  figure  below  for  an  illustration  of  the  mesh  pattern. ) 
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Figure    2.         N   =   5    >      M  =  4    . 
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As  in  the  one  dimensional  problem^  one  needs  to  express 
Wq  ,,    Vq  ,  as  a  linear  combination  of  the  values  at  neigh- 
boring points.   From  the  definitions  of  our  mesh  functions, 
we  have 

V     =   ¥ 
0,k      N,kM 

(3.3)  k  =  0,+l,+2, ...  . 

0,kM     l,k 

The  above  exhibits  the  overlapping  identification  of  some 
mesh  points  near  the  interface.   At  each  time  step  we  still 

have  to  prescribe  W^  rm+I  ^'^^  ^   "  0,+l,+2,...;  i  =  l,2, , 

M-1.   We  propose  to  do  this  in  terms  of  interpolation 
along  the  line  x  =  0.   In  general  one  could  prescribe 
such  points  as  a  linear  combination  of  the  V   .  values , 
say, 

+  0(Ay|+t+l) 

for  s,t  natural  numbers.   We  have  used  such  an  interpola- 
tion in  our  calculations  for  the  case  s  =  1,  t  =  2. 
Using  Taylor's  theorem  on  (3.4)  it  is  easy  to  show  that 
one  obtains  for  this  case 

_  _  ap(P+l)        ^       p(p+l)(a+l) 

(3.5) 

_   a(a+l)(p+l)           __    aP(a+l) 
"i,0  -       2  ^±,1 6 

where  a  =  i /M  and  p  =  1-a,  and  1  Jl  i  ^  M-1. 
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We  consider  a  mesh  refinement  problem  for  the 
Lax-Wendroff  [7]  9-point  difference  approximation  to  (3.I) 


(3.6) 


a 


t> 


)  v^ 


+1 


1    1    /i\ 
^1  fcl  ^'^    'j^,k^ 


a= 


'  «^^i^fefe^if^"5«,..,  J 


jA 


n  =  0,1,2,  .  .. 

J  =  -L,  1^)  ^,  .  .  . 

k  =  0, +l,+2j  . 


where  for  i  =  1,2, 


(3.7) 


^0,0 

c(i) 
^1,0 

c(i) 
^0,1 


2       2 

1  -  bf  .  -  bo  • 
1,1     2,1 

b^  .  +b,  . 


^2,i+^2,i 


^-1,0 


2 

b,  .  b, 
1 


c 


(i) 
0,-1 


\ 

i 

"i, 

1 

2 

^l 

i  ~ 

"2, 

i 

and  C_^/^  -  ^^\  ^'^    (jXk)  for  j  =  -l,lj  k=-l,l. 


where 


(3.8: 


l.i 


a  At 


Ax, 


and 


b 


2,1 


b  At 


Ay. 


To  treat  this  case  by  the  Krelss-Osher  results,  we  shall 
formulate  our  mesh  refinement  problem  as  a  system  of 
difference  equations  for  an  initial-boundary  value  problem. 
Because  of  the  more  complicated  nature  of  the  boiondary 
condition  in  two  dimensions,  we  found  it  necessary  to 
introduce  the  following  vector  of  dimension  M+1, 


-^0- 


(3.9) 


2"^  V 

JA 


n 

J  A 
n 

n 

j .  kM+1 


\ 


j,kM+M-l 


j  =  0,1,2,.. 
k  =  0,+l,+2, 


y 


¥e  now  consider  the  following  system  of  difference 
equations  equivalent  to  (3.6)a,b 


(3.10) 


,1     1    1 


A 


7^. 


j,k   ^^-^  1^-^   a,p  j-Ki,k+p 


for  n  =  0,1,2,...;  j  =  1,2,3,...;  k  =  0,+l,+2,  ... 
where  the  matrices    A   „   are  diagonal  and  given  by 


/ 


(3.11) 


A 


a,p 


.(1) 


\ 


.(2) 


v 


.(2) 


J 


The  above  approximation  (3.9)-(3.1l)  is  consistent 
with  the  following  system  of  partial  differential  equations. 
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(-  \ 


(5.12) 


w 


0 


W- 


1 


f~^^\ 


^ap^   0 
0    a 


V 


0 


0 


0 
a 


w 


0 
1 


X 


+ 


/bP2   0 
0    b 


\0 


0 


\ 


0 
Ob/ 


("    \ 


'0 


W- 


V^M-l/ 


y 


(where  e.g.  W^  ^^_^^  approximates  vj^(x,y,t)  =  u(x;,y+iAy-^^  t ) 
for  x  =  jAx^,   y  =  kMAy^  =  k  Ay^,   t  =  n  At). 

Formulating  (3.3)-(3.5)  as  a  matrix  boundary 
condition  we  obtain 


Vo,k 


(3.15) 


\ 


/ 


n 


W 


0,kM 


^0,kM+M-l 


\ 


0 


a 


1,0 


0 
0 
0 


a 


V 


M-1,0 


0 


o\ 

0 
0 


0 


^n 


a,k 
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/ 


+ 


^ 


0 
0 

a 


0 
0 


1,1 


M-1,1 


o 


V 

/ 


r 


\ 


0 
0 


1   0 


0 
0 


7^  + 

"^l.k+l 


O 


+ 


\ 


0 
0 
a 

a 


0 
0 


1,-1 


M-1,1 


V 


o 


,n 


'l,k-l 


( 

0 
0 

a 

a 


^n 


'N,k 


J 


\ 


) 


0 
0 


1,-2 


M-1,-2 


O 


V 


,n 


'l,k-2 


Theorem  V.   If  the  difference  schemes  (5.6)ajb  are  stable 
for  the  pure  Cauchy  problem  then  the  mesh  refinement 
scheme  (3.10)-{3.13)  is  a  stable  approximation.   (Note 
that  (3.6.)a,b  are  stable  for  the  pure  Cauchy  problem  if 

^1  1  -  ^/^  ^^'^   ^2  1  -  ^^'    ^^^    \1\') 

Proof.   To  apply  the  Kreiss-Osher  theorem^  we  must  find 

all  solutions  to  our  scheme  of  the  form 


(3.14)  z5^^=A"y^*j  ; 


y 


Ql 


1  ; 


Z 


<t>. 
-J 


<  00  . 


As  in  the  one  dimensional  case  of  constant  coefficients, 
we  consider  a  normal  mode  solution  to  (3.6)a 

=  1  . 


V.  ,  =  constant  A   i^  y^  ^ 


01 
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Then  T  satisfies  the  following  quadratic  equation^ 
■^,1"*-  "1,1  .  "1,1  "2,1       "1, 


(^L^Vi ,  ^.i>i , 


0 


"l.l  "2.1)^2 


0 


2 


(3.15)   +  ((1-^1^1-^2,1-^  +  (  ^^^2  ^'^)yo+( 


^2,1-^2,1 


2y 


0 


))^ 


^1  1  -^1  1    ^1  1  ^2  1 
+      2  5     ^0  ^ 


^1,1  ^2,1 


"W 


0 


0 


By  the  well  known  dissipativity  of  the  Lax-Wendroff 
scheme  [7],  an  application  of  the  same  argiiment  as  in 
Lemma  1  shows  that  for  [a]  _>  1,  A  ,/  1   that  the  x  roots 
are  separated  by  the  unit  circle.   Furthermore,  note  that 
the  case  y„  s  l  is  exactly  the  one  dimensional  case  treated 
in  Theorem  I.   The  above  proves  separation  of  the  modes 
for  T   ,x        the  respective  modes  of  equations  (3.6)a,b. 
This  allows  us  to  verify  stability  of  our  system  by  checking 
conditions  (a-c)  of  (2.10)  for  solutions  of  the  form 


(5.16)  Zj^^ 


^n  k 

A  yo 


Pii 


n  =  0,1,2,  ..  . 
j  =  0,1,2,... 
k  =  0,+l,+2,  .  .. 


where  the  p.  are  constants.   Substitution  of  (5.I6) 
into  our  boundary  condition  gives 
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4   -  R  T^ 


-2  ,       -1, 


(3.i7)<!p2  =  Po^i(«i,o+°^i,i^o  +  «i,-2yo  +«i,-iyo  ) 


Pm  =  Pn^i  (' 


From  the  above  equations  we  observe  that  a  non-trivial 


\ 


solution  will  exist  only  if  p„  and  p..  are  4   0.   This 
allows  us  to  eliminate  p^  and  p.,  in  the  first  two  equa- 
tions of  (3.17)  and  obtain 

T  tN  -  1 
1^2  -  ^  • 

As  in  Theorem  III  this  proves  stability,  since  the  diagonal 

elements  of  (3. 12)  -  ap-,  and  a  are  of  opposite  sign.  Q.E.D. 

Remarks . 

1)  The  above  argument  can  be  used  to  prove 
stability  for  a  general  boimdary  condition  of  the  type 
(3.5) •  For  any  such  condition  would  always  give  a  system 
of  equations  similar  to  (3.I7)  having  p^  and  p.,  as  the 
only  factors  on  the  right  hand  side. 

2)  As  in  Chapter  I^e,  we  note  that  it  is  possible 
to  extend  Theorem  5  to  the  matrix  equation 

u,  =Au  +Bu   !    A,Bj  diagonal  matrices, 
t      X      y      ' 

3)  The  analogous  two  dimensional  formulations  of 
Theorems  11,  IV  can  be  proven  by  the  above  method. 
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CHAPTER  IV.   COMPUTATIONAL  RESULTS 

Calculations  were  performed  using  the  mesh  refine- 
ment schemes  of  Theorems  I-III,  and  V.   The  primary  purpose 
was  to  find  out  the  size  of  the  inaccuracies  which  are 
introduced  by  the  different  interface  conditions  discussed 
above.   Due  to  the  hyperbolic  nature  of  our  problem,  we 
cannot  in  general  expect  to  achieve  higher  order  accuracy 
in  the  refined  regions.   This  is  so,  because  the  error 
from  the  unrefined  region  can  be  expected  to  propagate 
into  the  refined  region.   In  fact,  our  calculations  were 
performed  for  one  period  of  a  periodic  problem.   Because 
of  this,  all  regions  will  eventually  exhibit  the  same 
order  of  accuracy. 

The  mesh  refinement  technique  might  be  a  practical 
method  for  handling  problems  where  large  gradients  in  the 
solution  are  caused  in  a  confined  region.   In  such  cases, 
if  new  inaccuracies  do  not  arise  at  the  interface  of  the 
mesh  patterns,  the  mesh  refinement  scheme  enables  us  to 
efficiently  make  a  local  adjustment  to  obtain  higher 
accuracy. 

a .   One  dimensional  periodic  case. 

For  our  first  series  of  calculations,  we  solved  the 
initial  value  problem  of  Example  I  with  a  s  1  and  the 
solution 
(4.1)  u(x,t)  =  sin  hv{x   +t)  . 
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We  used  the  periodicity  of  our  solution  to  restrict  our 
calculations  to  the  interval  [0,1].   We  refined  the  region 


(^.2) 


ij^  =  {x  :  J  1  X  ^  J  3  . 


The  mesh  sizes  were  such  that 
Ax. 


(^•5)        p  =  s^  =  ro^ 


1 


^^2  "   150  ' 


At 


1 


1750  • 


n      n 
Let  V.  and  w.  denote  the  unrefined  and  refined  schemes, 

respectively.   We  define  (see  Figure  3  below) 


(^.4)     v^i  -   ™1  > 


n      n 

^101  =  ^501  ' 


and  by  periodicity  we  always  have 


(4.5) 


n    n 

^1  =  ^151 


V-, 


— » — \ 


^101=  ^501 


►I- 


'151 
- — I 


x=0 


-h 


2 
^=5 


x=l 


Figure  3.   Refined  Middle  Third  of  Unit  Interval. 

Of  course,  we  do  not  compute  v.    for  i=  52,53, ... ,100, 
but  use  the  values  of  w . 

We  describe  the  accuracy  of  the  mesh  refinement 
scheme  by  computing  several  quantities.   Let  e(v)  and  e(w) 
denote  the  maximum  error  of  our  scheme  from  the  solution 
in  the  \inrefined  region  and  the  interior  of  I„,  respectively. 
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The  size  of  the  error  near  the  interfaces  depends  on  the 
direction  of  characteristics.   In  our  case,  the  error 
from  the  imrefined  region  advances  into  I  at  x  =  2/5. 
To  the  right  of  x  =  1/5  our  calculations  show  an  abruptly 
improved  accuracy  until  the  lonrefined  solution  has 
propagated  up  to  x  =  1/5  from  beyond  x  =  2/5.   This  happens 
after  approximately  5OO  time  steps. 

To  illustrate  what  the  error  in  I„  would  be  if 

K 

the  characteristics  always  emerged  from  I  ,  we  also  compute 
a  quantity  e  .  (for  £  =   lOO)  defined  by 

(4.6)  e^=  max     |u(x,t)  -  w^|  ;   t=nAt,  x=  ^  +{j-1)^k^. 

In  Table  1,  we  list  the  results  of  several  calcula- 
tions for  the  example  of  (4.1)  using  the  Lax-Wendroff 
scheme  with  refinement.   The  entries  in  rows  I-5  exhibit 
the  stability  and  relatively  improved  accuracy  in  I„  of  the 
mesh  refinement  scheme.   The  three  different  interface 
conditions  give  nearly  identical  accuracy. 

The  entries  in  rows  4  and  5  correspond  to  results 
using  no  refinement  with  Ax  =  Ax„.   For  row  4  we  used 
At  =  1A750-   Fo^  ^ow  5  we  took  our  time  step  as  At'=  lA75 
and  thus  iterated  only  40  times.   It  appears   that  one 
loses  accuracy  by  using  too  small  a  time  step  (aside  from 
having  to  do  more  work) .   This  result  is  in  line  with  a 
commonly  observed  fact  that  one  gets  better  results  when 
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Table   1 


a  =   constant 


At   =   IA750;         At'    =   IA75;         ^x^   =   IA50;         Ax-j^   =   IA5OO 


t,  time 


e(v) 


System  I 

System  II 

System  III 

No  Refinement 
(Ax^) 

Uneven  Time  (l) 
(III) 

Unstable  (ll) 

Non-Dissipative 


h 


At^ 


4X10 


-6 


400  At 
4oo  At 
4  00  At 
400  At 

400  At' 

40  At' 
40  At' 

4  00  At 

5000  At* 


3.332X10  ^ 
3.332x10"^ 


3.332XLO 
1.679x10 

8.909x10 

8.909XLO 
8.909x10 

2.461 

4.624XL0 


-3 
-3 


-4 


-4 


-4 


-4 


e(w) 


I.59IXLO  ^ 
1.656x10"^ 


1.629x10 


-3 


4.493XLO 
4.493x10 

3.190 

1.752x10' 


-4 


'100 


6.3x10 
6.3x10 
6.7x10 


-6 
-6 


-4 


-4 


6.3x10' 
6.3XLO" 


-6 


3.3x10  ^    i 
I 

1 
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the  stability  factor  bp  =  a  At/Ax„  is  close  to  one  in 

absolute  value.   In  the  mesh  refinement  problem,  a  small 

time  step  is  needed  for  stability  in  I  .   However,  using 

the  same  time  step  in  the  imrefined  region  will  reduce 

accuracy  there.   To  rectify  this  situation  we  also 

considered  using  different  time  steps  in  the  different 

regions.   In  I_  we  used 
a 

(4.7)  At  =  1A750 


as  above.   In  the  xinrefined  region,  we  used  a  time  step 
(4.8)  At'  =  ^  =  ^ 


P  "  175  • 

The  only  difficulty  in  using  this  approach  is  at  the 
interface.   We  describe  how  we  treated  the  interfaces 
at  X  =  1/5  for  Systems  I  and  III.   A  similar  fomiulation 
was  used  at  x  =  2/5. 

Each  cycle  of  computation  consists  of  first  advancing 
wj^"^-^  and  w^/P  '^^   for  i=l,2,  ...10.    Then  the  remaining 
is  advanced  in  the  usual  manner  until  all  points  are  at 
the  same  time  level.   For  System  I  we  used  the  following 
modification  of  (1.21) 


i4.9ja  w^     _  [T^  ^50 


^P  i  n      ^P  i  (^P  i+1)    r^  M 


where 


5i> 


(4.9)b         b^^.  =  ^i^   and   p.  =  ip 

for  i=  1,2,  ...,10. 

For  System  III  with  b„  .  as  above,  we  used 

Rows  6  and  7  of  Table  1  show  the  results  obtained.   Aside 
from  the  additional  labor  saved,  the  main  advantage  of  the 
uneven  time  step  approach  is  the  improved  accuracy  outside 
the  refined  region. 

In  a  remark  after  Theorem  II,  we  stated  that  one 
must  consider  the  domain  of  dependence  of  the  solution 
when  choosing  an  interface  condition.   To  illustrate  the 
instabilities  which  may  arise  when  Theorem  II  is  not 
observed,  we  calculated  the  scheme  of  row  2  as  above  with 
a  =  -1.   This  has  the  effect  of  reversing  the  direction  of 
characteristics.   The  interface  conditions  however  were  left 
unchanged.  In  row  8  the  errors  listed  clearly  indicate 
an  instability. 

As  a  further  example,  we  calculated  System  II  above 
with  a  different  underlying  difference  scheme  namely 

(4.11)   «-l  .  1^1  „n.^  .  „n  ^  ^  „„^^  _ 

Our  theory  does  not  apply  to  this  scheme,  since  it  is 
non-dissipative,  does  not  have  separation  of  the  modes, 
and  in  fact  is  stable  for  a  At/Ax^  <_   constant.   Despite 
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this,  the  results  of  our  calculations  in  row  9  indicate 
that  the  mesh  refinement  scheme  is  stable. 

In  Table  2  we  present  the  results  of  computations 
performed  for  Example  I  with 


(4.12)   a  =  a(x) 


The  analytic  solution  can  be  found  by  the  method  of 
characteristics  so  that 


1 

2 

2  +  s  in  Gvx 
4 

^  1  ^'  ik 

X     6     1^ 

1 
2 

2 

^    <_  X    _<    1 

u 


{x,t)    =      <l>(x)        where      '^{x)    =   sin  hirx 


(^.13)  5^ 

d^ 


and  t   = 


The  quantities  displayed  in  Table  2  are  defined 
exactly  in  the  same  way  as  in  Table  1.   We  also  computed 
with  several  time  steps  for  no  refinement  and  uneven  time 
steps.   We  note  that  in  this  case  it  was  not  immediately 
obvious  whether  the  uneven  time  steps  would  improve 
accuracy  in  the  unrefined  region. 
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Table   2 


a  -  a(x)    as   in    (4.1) 
At  =  IA250;        At'    -   IA25;        Ax^  -  IA50;        Ax-^  =  IA5OO 


1)   System  I 

1 

tf    time 

e(v) 

e(w) 

^100 

400  At 

7.136X10"^ 

1.397X10"^ 

l.OOOXLO"^ 

2)    System  II 

400  At 

7.134X10"^ 

2.764X10"^ 

l.OOOXLO"^ 

5)    System   III 

400  At 

9.993x10"^ 

1.401X10"^ 

1.000X10"^ 

4)   No  refinement 
(Axg) 

400  At 

l.l8oxio~^ 

5)    " 

40  At' 

2.358x10"^ 

6)    No  refinement 

27  At* 

2.901x10"^ 

At*   =  I5A25O 

7)    Uneven   time    (l) 

40  At' 
=   400At 

7.668x10"^ 

1.652x10"^ 

l.OOOXLO"^ 

8)   Uneven   time  (iTf) 

40  At' 

7.668x10"^ 

1.6llXL0"^ 

i.oooxLo"^; 

=   400At 

1 

9)    " 

27   At* 

4. 841X10"^  jl .  553X10"^ 

1. 086x1 0"^j 

=  405At 

J 

i 

1 

) 

1 
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t).   Two  dimensional  periodic  case. 

Calculations  were  performed  for  the  initial-value 
problem 

^t    =   -  ^x  +  ^y 
(4.14) 

u(x,y^O)  =   <t)(x,y)  =  sin  2fx  cos  2Try  . 

We  used  the  periodicity  of  the  solution 

u(x;,t)  =  sin  27r(x-t)  cos  27r(y +t) 

to  restrict  our  calculations  to  the  unit  square 
0  <  x,y   <_   1.   The  refined  region  was  taken  as  (see 
Figure  4  below ) 

(^!.15)         DR={(x,y):  |ix,y<|3. 

Using  the  notation  of  Chapter  III  we  took 

Ax„  =  Ay^  =  ^  ;    Ax.  =  Ay^  =   ^ 


-2  =  "^2  =  4:5  ^    "^1  =  "^1  -  225 
(^.16) 

Pi  =  Po  =   1-    I         ^t 


1  ~  ^2  ~  5  '  1750  • 

As  in  Table  1,  we  describe  the  accuracy  with  the 
analogous  quantities  e(v)  and  e(w),  the  maximum  errors  of 
our  scheme  in  the  unrefined  region  and  D^,  respectively. 
To  describe  the  accuracy  of  the  refined  scheme  away  from 
the  interfaces,  we  compute  the  error  in  the  interior  of  D„, 
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1,46 


""he.he 


« 
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« 

t         » 

.^l,76=''l6,36 

1^^76,76=^31,51 

°R 

• 

g 

• 

•  ••»«. 

•  •  «  -  • 

• 

*"'****"* 

^1^=^16,16 

"16, 

1=^51,16 

•             •           • 

•                ••              • 

•             •           • 

*       • 

»            •          • 

•           *          * 

V 


1,1 


'46,1 


Figure  4.      Unit   square,    0  j<  x,y  j<  1. 
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(4.17)   e.  ,   =    max    |u(x,t)-w^   |  ;    y  =  l+(k-l)Ay 

t  =  n  At 


X  =  J  +(j-l)Ax^ 


In  Table  3  we  list  the  results  of  our  calculations 
using  the  Lax-Wendroff  scheme  (3-7)  with  our  mesh  refine- 
ment scheme  of  Theorem  V.   We  compare  the  above  with 
ordinary  Lax-Wendroff  with  three  choices  of  time  steps. 
We  then  compute  with  these  time  steps  using  uneven  time 
steps . 

The  procedure  for  computing  with  uneven  time  steps 
is  similar  to  the  one  dimensional  case.   We  first  compute 
the  interface  by  using  the  two  dimensional  analog  of  (4.10) 
at  all  v^   on  the  interface.   The  only  difficulty  in  using 
(5.5)  is  for  points  near  the  corners.   At  uneven  time  levels 
when  interpolating  for  the  values  w.  ,  where  2  ^  j  <_  5  or 
72  <  k  <  75,  formula  (3.5)  calls  for  values  of  v .  ,   which 
are  not  available.   In  such  cases  we  substituted  a  3  point 
interpolation  scheme  again  of  the  type  (3.4). 

Our  computations  show  that  the  interface  condition 
(3.5)  does  not  introduce  any  additional  errors.   When 
uneven  time  steps  are  used,  the  error  in  the  unrefined 
region  compares  with  the  error  for  ordinary  no  refinement. 
While  for  all  entries,  ^±Yxt      remains  the  same.   Thus,  one 
should  try  to  separately  adjust  the  time  steps  so  as  to 
minimize  the  error  in  each  region,  and  then  match  the 
schemes  along  the  interface. 
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Table   5 


At    =   IA750;         Ax^  =  Ay^   =    IA5;         Ax^  =  Ay-^_   =   1/225 


1)  Ordinary 
refinement 

2)  No  refine- 
ment 
(Ax^jAyg) 

a)  At  =  IA75O 

b)  At'=   5A75O 

c)  At*=  1A75 

3 )  Uneven  time 


tj    time 


100  At 


100  At 


20  At' 


10  At 


20At'=100At 


lOAt   =100At 


e(v) 


I.I82KLO 


-3 


l.l65>^0"^ 
1.104x10"^ 
8.895x10"''^ 

1.104x10"^ 
8.895x10' 


e(w) 


1.143x10"^ 


1.033x10 
7.822x10 


-3 

-4 


'int. 


4.36x10 


-5 


4.36x10"^ 


4.36x10 


-5 
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c.   Conclusions . 

In  this  chapter,  we  presented  numerical  results 
which  seem  to  Indicate  the  stability  of  the  mesh  refinement 
technique  for  several  cases  not  analytically  treated  above. 

Our  main  observation  was  the  Improved  accuracy 
gained  by  using  uneven  time  steps  with  the  mesh  refinement 
scheme.   A  stability  proof  of  suQh  a  scheme  Involves  an 
Imposing  amount  of  algebra.   In  the  future  we  hope  to  prove 
stability  only  for  a  simple  case. 

¥e  shall  also  try  to  extend  these  techniques  to 
solving  parabolic  equations. 
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